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Abstract 

This paper introduces a new Windowed Green Function (WGF) method for the numerical 
integral-equation solution of problems of electromagnetic scattering by obstacles in presence of 
dielectric or conducting half-planes. The WGF method, which is based on use of smooth window¬ 
ing functions and integral kernels that can be expressed directly in terms of the free-space Green 
function, does not require evaluation of expensive Sommerfeld integrals. The proposed approach 
is fast, accurate, flexible and easy to implement. In particular, straightforward modifications of 
existing (accelerated or unaccelerated) solvers suffice to incorporate the WGF capability. The 
mathematical basis of the method is simple: the method relies on a certain integral equation 
posed on the union of the boundary of the obstacle and a small flat section of the interface 
between the penetrable media. Numerical experiments demonstrate that both the near- and 
far-held errors resulting from the proposed approach decrease faster than any negative power 
of the window size. In the examples considered in this paper the proposed method is up to 
thousands of times faster, for a given accuracy, than a corresponding method based on the 
layer-Green-function. 


1 Introduction 

The solution of problems of scattering by obstacles or defects in presence of planar layered dielectric 
or conducting media has typically required use of Sommerfeld integrals and associated layer Green 
functions—which automatically enforce the relevant transmission conditions on the unbounded flat 
surfaces and thus reduce the scattering problems to integral equations on the obstacles and/or 
defects. As is well known, however, the numerical evaluation of layer Green functions and their 
derivatives, which amounts to computation of certain challenging Fourier integrals [8, 20], are ex¬ 
tremely expensive and give rise to a significant bottleneck in layer-media simulations (see e.g. [6] 
for details). This paper presents a novel integral-equation approach for problems involving layered 
media. The new approach, which is based on use of certain “windowing” functions and considera¬ 
tions associated with the method of stationary phase, does not require use of expensive Sommerfeld 
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integrals. Numerical experiments demonstrate that both the near- and far-held errors resulting 
from the proposed approach decrease faster than any negative power of the window size. 

A variety of methods have been provided for the solution of problems of scattering by obstacles 
in presence of layered media. Amongst the most effective such approaches we mention 1) Methods 
which evaluate Sommerfeld integrals on the basis of path-integration in the complex plane [17, 7, 
6, 18] (such approaches require numerical evaluation of integrals of functions that oscillate, grow 
exponentially in a bounded section of the integration path and, depending on the relative position 
of the source and observation points to the interface between the two media, may decay slowly at 
inhnity); 2) The complex images method reviewed in [1] (a discussion indicating certain instabilities 
and inefficiencies in this method is presented in [7, section 5.5]); and 3) The steepest descent 
method [9, 10] which, provided the steepest descent path is known, reduces the Sommerfeld integral 
to an integral of an exponentially decaying function (unfortunately, however, the determination of 
steepest descent paths for each observation point can be challenging and expensive). As is well 
known, in any case, all of these methods entail significant computational costs [6]. 

The approach proposed in this paper bears similarities with certain “finite-section” methods 
in the field of rough-surface scattering. These methods utilize approximations based on truncated 
portions of a given unbounded rough surface [14, 22, 19] and, in some cases, they incorporate 
a “taper” [22, 21, 15] to eliminate artificial reflections from the edges of the finite sections. In 
fact the smooth taper function utilized in [15] (Figure 2 in that reference) resembles the smooth 
windowing function we use (Figure 2 below and reference [3]). But as indicated in comments 
provided in section 2 below in regards to certain slow-rise windowing functions, essential differences 
exist between the finite-section approaches and the methods proposed in this paper. In particular, 
with exception of the slow-rise windowing function method [3, 16], none of the previous tapered 
rough surface algorithms has demonstrated high-order convergence as the width of the finite sections 
tend to infinity. 

In section 4 the proposed WGF method is compared against the high-order integral equation 
method recently introduced in [18], which is based on the accurate and efficient evaluation of the 
Sommerfeld integrals. In the examples considered in that section the proposed method is up to 
thousands of times faster, for a given accuracy, than a corresponding method based on the layer- 
Green-function. A much larger improvement in the computational cost is expected for problems 
of electromagnetic scattering by defects and obstacles in multi-layer structures in two- and three- 
dimensional spaces, which will be addressed in future contributions. 

The proposed methodology is presented in sections 2 and 3. A variety of numerical results 
presented in sections 2 and 4 demonstrate the accuracy and speed of the proposed approach. 

2 Windowed Green Function Method 

We consider two-dimensional TE and TM polarized dielectric transmission problems. As is well 
known, the z components u = and u = Hz of the total electric and magnetic fields satisfy the 
Helmholtz equation Au+kjU = 0 in Qj, j = 1,2 (see Figure 1), where, letting w > 0, Sj > 0, /xq > 0, 
and aj > 0 denote the angular frequency, the electric permittivity, the magnetic permeability of 
vacuum, and the electrical conductivity, the wavenumber kj is defined by k'j = + iaj/uj)fio, 

j = 1,2. In either case the total field is given by 

f in Hi, , . 

[ U 2 m U 2 , 
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where denoting by a G (—7r,0) the incidence angle measured from the horizontal (see Figure 1), 
^mc(^) _ g*A:i(xi coso+a; 2 sma)^ denote the incident plane-wave and the reflected and 

transmitted waves, respectively. As is known (see e.g. [11]), the scattered and transmitted fields ui 
and U 2 admit the representations 


ui = Vi [(f \ — 5i [V'j in 111 , (2a) 

U2 = -^>2 M + ^2 [V’j in fl2, (2b) 

in terms of the total field = u\y and its normal derivative ■0 = on F, where letting Gj{x, y) = 
— 2/|)/4, j = 1,2 denote the free-space Green function for the Helmholtz equation with 
wavenumber kj, the single- and double-layer potentials in equation (2) are defined by 


•SjblKx) = j Gj{x,y)rj{y)dsy, and 


Vj[y]{x) = 


'r duy 


{x,y)y{y)dsy, 


(3) 


respectively. By evaluating the fields (2) and their normal derivatives on F and using the transmis¬ 
sion conditions 

jjjj, du2 dui 


U 2 — Ul = u 


dn dn dn 


(with V = 1 and v = ei /£2 in TE- and TM-polarizations respectively) we obtain the second-kind 
system of integral equations [12] 

E0-bT0 = 0'“ on F (4) 

for the surface currents 0, where 


E = 


1 0 

0 4^ 
^ 2 


and where 


T = 


u|r 

I 9u I 

L 9n 


D2 - Di -uS2 + 5i 

N2-N1 -VK2 + K1 


dn n 




/ dG 

^{x,y)y{y)dsy 


(5) 


is defined in terms of the boundary integral operators defined by the expressions Sj[rj\{x) and 
Dj[ri\{x) as well as 


for a; G F and for j = 1, 2. 

Instead of solving the problem on the entire infinite plane a locally windowed problem could 
be used in an attempt to obtain local currents over relevant portions of the geometry. To pursue 
this idea we introduce a smooth windowing function wa (which is depicted in Figure 2) which 
is non-zero in an interval of length 2A, and which has a slow rise: wa{xi) = f{xi/A) for some 
fixed window function /. (Note that, with such a definition, wa rises from zero to one in a region 
of length proportional to A; see [3, 16]. As demonstrated in those references, the slow rise of 
the window function is essential to ensure fast convergence of the approximation.) For notational 
simplicity, the subindex A will be dropped in what follows, and we will thus write rc(xi) instead of 
wa{xi). The parts of the boundary F where rc(xi) 0 and w{xi) = 1 — w{xi) 0, further, will be 
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Figure 1: Description of the problem under consideration: scattering by a defect in a dielectric or 
conducting plane. F denotes the interface between the two media while IT denotes the interface 
between the upper- and lower-half planes. 


denoted by F^ and F^i, respectively. The width 2A > 0 of the support of the window function w is 
selected in such a way that tc(a:i) vanishes on any corrugations that exist on the surface F, as well 
as on any additional obstacles that may exist above and/or below F. (For notational simplicity 
our derivations are presented for cases for which the corrugations on the surface F are the only 
departures from planarity, but, as demonstrated by Figure 12, our algorithms are also applicable 
in cases in which additional scatterers exist.) 



Figure 2: Window function w = wa and the windowed sections F^i and F^ of the unbounded 
curve F. 

Utilizing the windowing function w and letting W = w ■ I, where I is the 2x2 identity matrix, 
we consider the preliminary approximate equation 

E(j)*+ TW(j)* = on Ta (6) 

(where the new unknown cf)* is defined on F^), and, in order to assess the errors inherent in this 
approximation, the form 

Ecj) + TWcj) = />'“ - T{I - W)(j) on Ta (7) 

of the exact equation (4). Using integration-by-parts and employing the method of stationary-phase, 
it follows [5] that the term T{I — W)4> is super-algebraically small (i.e., smaller than Cp{kA)~P for 
any positive integer p as kA —)• oo, where Cp is a p-dependent constant) in the region {tc = 1}, and, 
thus, as shown in [5], that the solution cjA of (6) is a highly accurate approximation of cj) throughout 
the center region {w = 1} of the surface F^ provided A is large enough. However, it is easy to 
see that, to correctly take into account fields reflected from the planar portions of the surface, the 
needed window sizes may be very large—especially so for incidence angles approaching grazing. 

To demonstrate this fact we use equation (6) to approximate the solution of the TE problem 
of scattering of a plane-wave by a semi-circular bump of radius a = 1 placed directly on top of a 
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planar dielectric surface. The problem was discretized using a graded mesh over the surface of the 
bump and on the windowed portion of the planar interface, on the basis of a direct generalization 
of the Nystrom method presented in [13] with p = 3. For this example the wavenumbers ki and k 2 
in the regions above and below the plane were set to dvr and Svr, respectively, and approximately 
20 points per unit length of the surface of the bump and the surrounding were used. 

As shown in Figure 3, the naive windowing approach embodied in (6) requires large regions of 
the planar interface to be discretized as the incidence angle decreases. For accurate calculations at 
even moderate angles, a large number of wavelengths must be present in the window region, well 
beyond the extent of the non-planar local geometry. 




Figure 3: Errors in the integral densities resulting from numerical solution of (6) by means of 
a naive implementation of the WGF method for a semi-circular bump-shaped defect, for various 
window sizes and angles of incidence. Left: log-log scale. Right: semi-log scale. Clearly, the window 
size required by the naive method to produce a given accuracy increases dramatically as the angle 
of incidence approaches grazing. 

In order to provide an insight into the source of the errors displayed in Figure 3 we present 
Figure 4. Figure 4(a) presents rays incident on the left planar region as well as their reflection and 
transmission. Clearly, in view of the incidence angle considered these reflected fields subsequently 
illuminate the defect. The blue rays, for example, represent the reflections that are correctly taken 
into account in the solution of equation (6) (since they lie within the windowed region), while the 
red arrows represent reflections that are neglected. Figure 4(b), on the other hand, represents 
reflections by the defect. The color-code in the left figure carries over to the right figure: the blue 
(resp. red) rays in Figure 4(b) represent the fields scattered by the defect which arise from the blue 
(resp. red) arrows in Figure 4(a). We remark that the scattering of the field represented by the red 
arrows is not taken into account by (6), which gives rise to the errors observed in Figure 3. We also 
note that the relatively fast convergence demonstrated by the blue curves in Figure 3 is explained 
by the fact that for near normal incidence {a ~ —7r/2) there is not much “red field” interacting 
with the defect. In contrast, for incidence near grazing (a ~ 0), “red fields” from regions far away 
from the windowed area do interact with the defect. This explains the poor convergence properties 
demonstrated by the green and red curves in Figure 3: the fields neglected in the naive approach 
give rise to important contributions as a decreases. 

To address this difficulty we consider again the exact integral equation (7) and we substitute the 
unknown density </> on the right-hand side of this equation by the corresponding (known) density (j)^ 
associated with the problems of scattering and transmission of a plane-wave by a perfectly flat 
infinite plane. Since a super algebraically small portion of the field reflected by the windowed region 
reflects back into the windowed region upon reflection from the plane outside the windowed region. 
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we conclude that the error arising from the substitution of (j) by results in superalgebraically 
small errors in equation (7) throughout the region {w = 1}. We thus obtain the approximate 
equation 

+ TWcT = - T{I - on Fa, (8) 

whose solution 0"' is a superalgebraically close approximation of the exact solution (p throughout 
the region {w = !}• In order to evaluate the term T{I — W)4>^ we note that since {I — W)(j)^ is 
zero everywhere deviates from the planar boundary 11 = {(xi,X 2 ) € : X 2 = 0} (depicted in 

Figure 1), we have 

T{I - W)pf = Tn{I - W)p^, 
where letting the layer potentials 5? and be given by 


•Sf[v]{x)= / Gj{x,y)rj{y)dsy, and 
J n 

^-^{x,y)y{y)dsy, 

the operator Tn is defined as 

^ r -uS^ + SY ■ 

^ [ A^2" - iV" 


(9) 


in terms of the boundary integral operators defined by the expressions SY[ri]{x) and DY[ri]{x) as 
well as ^ 

= J^^^ix,y)y{y)dsy 

for a; E r and for j = 1, 2. Thus equation (8) becomes 

+ TWP^ = - Tnp^ + TuWpf on F^. (10) 


Clearly the expression TjiW can be evaluated by means of integration on the bounded region 
n n {(xi,X2) E : w{xi) 7 ^ 0}, and the expression can be computed in closed form: 


Tupf = 


/ 


L.nc 

T 



a a , 

on 



„,inc 

J - {l + u)uf/2y 

y. 


dn 



on F \ n, 


on F n n, 


( 11 ) 


where is the total field resulting from the solution of the problem of scattering by the flat 
dielectric plane with boundary FI [8, Chapter 2], 
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From the discussion above we see that, on the set {re = 1}, the (superalgebraically high) 
accuracy of the solution (j)'^ of (10) does not deteriorate as the incidence angle a tends to zero. 
As shown in section 3 below, further, the solution 0"' can be used to produce the total field u 
everywhere in space as well as the associated far field pattern. To conclude this section, in Figure 5 
we demonstrate the fast and angle-independent convergence of 0"' to cj): clearly the value of A 
required to obtain an accurate approximation of the exact solution has been reduced substantially 
and the errors are uniformly small as the incidence angle decreases to zero. 



A/\ 



A/\ 


Figure 5: Errors in the integral densities on the surface of the defect resulting from numerical 
solution of (10), for a semi-circular bump-shaped defect, and for various window sizes and angles 
of incidence—including extremely shallow incidences. Left: log-log scale. Right: semi-log scale. 
Clearly, this version of the WGF method computes integral densities with super-algebraically high 
accuracy uniformly for all angles of incidence (cf. Figure 3). 


3 Field evaluation 


An analysis similar to the one presented in section 2 for the density cj)^ = shows that 

substitution of (/> = by + (1 “ w)tf ^+ (1 “ w)'ij)fY ™ (2) produces the fields ui 

and U 2 with superalgebraically high accuracy in a neighborhood of the region {w = 1} in and, 
in particular, on a closed disc D such as the one depicted in Figure 7. After some manipulations 
similar to those presented in the derivation of (11) above, the resulting formula can be re-expressed 
into a formula for the total field in terms of surface potentials defined on both F and Ft, namely 


u{x) = Vi [vinp^] (x) — Si (x) — 

u^{x), X G {X 2 > 0}, 


Wif’ 


(x) + 5]^ (x) 


(12a) 


+ 


0, ® G {x2 < 0} 


for a; G Hi, and 


u{x) = -V2 [wip'^] [x) + S2 {x) + vf 

0, a; G {x2 > 0}, 


Wif' 


(x) - 


vw 


Ipf ( 


(X) 


(12b) 




I'^ix), X G {X2 < 0} 


for X G ^^2- 
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Figure 6 compares the total field obtained by means of the WGF method and the layer-Green- 
function method [18] for the solution of the problem of scattering of a plane-wave by a semi-circular 
bump of radius a = 1 in TE-polarization for wavenumbers ki = 10 and /c 2 = 15 for a = —7r/2 
and OL = —tt/G incidences. The WGF solution, in particular, was obtained from the solution of 
the integral equation (10) followed by evaluation of field values on the basis of (12). Figures 6c 
and 6f, which display the absolute value of the difference of the total fields computed using the 
WGF method and the layer-Green-function method on a bounded portion of the strip {tc = 1} 
demonstrate the accuracy of the computed solutions in the near field. 





(a) WGF method, 
s 
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-6 

-8 



(b) LGF method. 



-2 0 2 


(c) Difference. 

in-l 



-2 0 2 


(d) WGF method. 


(e) LGF method. 


(f) Difference. 


Figure 6: Real part of the total fields produced by the WGF method (first column) and the 
layer-Green-function method [18] (second column), and absolute value of the difference of the fields 
computed using the WGF method and layer-Green-function method (third column) for the problem 
of scattering of plane-wave by a semi-circular bump for a = —7r/2 (first row) and a = —vr/G (second 
row) incidences. The width of the support of the selected window function is 2A = IGA ~ 10.053 
in all these calculations. The black lines represent the domains of the respective integral equation 
formulations. 


As may be expected, however, formulae (12) do not generally provide an accurate approximation 

























































of either far fields or near fields outside a neighborhood of Fyi. In order to tackle this difficulty 
we consider the boundary S of the disc D mentioned above and depicted in Figure 7: S encloses 
the portion of F that differs from the flat interface FI and, as indicated above, it lies within a fixed 
region within which super algebraic convergence of the fields ui and U 2 takes place. Application of 
the Green identities, integrating over the region exterior to S and utilizing the layer Green function 
leads to the following integral representation of scattered field = u — u^: 

r (pjr..s 'I 
u"{x) = - G\{x,y) — {y) | dsy (13) 

outside the region enclosed by 5, where G 2 denotes the layer Green function for the Helmholtz 
equation with wavenumbers ki in {x 2 > 0} and k 2 in {x 2 < 0} that satisfies homogeneous trans¬ 
mission conditions on the flat interface H (see Appendix A). Note that the scattered field and its 
normal derivative on S can be computed directly utilizing (12) since by construction S lies inside 
the region where (12) provides an accurate approximation of the total field u. 



Figure 7: Surface S utilized in (13). 


The far-held pattern ttoo 
u^x) = 


(x), which is related to the scattered held by the asymptotic formula 

— ^Uoo{x) + r = |£c| ^ 00 , x = —, 

y/r \x\ 


can be obtained from (13) in a straightforward manner by replacing G 2 by its asymptotic expansion 
as I a; I — >■ 00 . The hrst order term of the asymptotic expansion of the Sommerfeld integrals <I>i 
and < 1>2 (equation 20) in a given direction x = (cos a, sin a), 0 < a < n can be obtained by the 
method of steepest descent by taking into account the contribution of the saddle point [9] (branch 
point singularities and poles do not contribute to the hrst term of the asymptotic expansion of the 
two-layer Green function). Substitution of the result in equation (13) gives rise to the expression 


'^ooix) — / 

Js 


'x, y)u^{y) - H{x, y) — {y) \ dsy 


dn 


dn 


for the far held Uoo{x), where 


H {x,y) =■ 


v{k2 — k\) ^-ik]_x-y+i’K/ a 


V2vrA:i(l -F v) (m + Vi) ivi + 
for y E { 2/2 > 0} and 




VSvrfci 




1 — e 


l + v 


^-ikix-y+i'n'/A 


x/Svr/ci 


(14) 


(15a) 


H {x,y) 


vki sin(a - (3) 
y/27rki rji + vri2 


(15b) 
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for y G {7/2 < 0}, where S = x/\x\ = (cosa, — sina), y = |7/|(cos/3,sin/?), rji = 7 i(/ci cos(q; — /?)) 
and 772 = 72(^1 cos(q; — /3)) (see Appendix A for the definition of 71 and 72 ). Thus, unlike the 
layer Green function G\ itself, the far field associated with G\ can be computed inexpensively by 
means of the explicit expressions (15). Figure 8 provides a comparison of the far-held patterns 
computed using the layer-Green-function method and the WGF method proposed in this paper for 
the example problem considered above in the present section 3. 


9° 0.5 

—WGFM 

90 1.5 

WGFM 

120 60 

. LGFM 

120 60 

. LGFM 



Figure 8: Far-held patterns obtained using the layer-Green-function method [18] (red dotted curve) 
and the WGF method (continuous blue line) for the solution of the problem of scattering considered 
in this section at incidences a = —7r/2 (left) and a = —vr/b (right) . 

In view of this discussion, equations (12) and (13) can be used to accurately and efficiently 
evaluate near-helds and far-helds, respectively. These are typically the quantities of interest in 
scattering simulations involving layered media. The evaluation of the helds in an intermediate 
region, such as a domain outside the neighborhood of F^ where (12) yields an accurate approxi¬ 
mation, can also be approximated efficiently on the basis of equation (13). Indeed, in such cases, 
for which source points y lie on S and observation points x are at a certain distance away from 
S', the Sommerfeld integrals (20) and (22) (which contain highly oscillatory and/or exponentially 
decaying integrands) can be obtained by means of asymptotic numerical methods [2, 4] based on 
localization around critical points [9, 18]. 

4 Numerical Experiments 

This section illustrates the proposed methodology with a variety of numerical results concerning 
dielectric and conducting media, including relevant efficiency and accuracy studies. 

In our first example we consider once again the configuration associated with Figure 5 (i.e. the 
problem of scattering by a semi-circular bump defect on a dielectric plane in TE-polarization). Here 
we compare the computing times required to create the systems of equations (which is the operation 
that dominates the computing time in all the examples considered) that stem from the discretization 
of the relevant integral equations by means of the WGF method (10) and the layer-Green-function 
method [18, Eq. 7]. Eigure 9 displays the computing times for various wavenumbers fci and ^2 = 2/ci 
for each method. The discretization density was held proportional to ki to properly resolve the 
oscillatory character of the integrands and the same discretization was used for both methods on 
the bump, allowing for a point by point comparison of the solutions. In all these examples the WGF 
method was optimized to produce a maximum error of approximately 5 x 10“^ in the computation 
of the density (j)^ on the surface of the bump. Similarly, the key parameters in the implementation 
of layer-Green-function method (including the parameters associated to the numerical evaluation 
of the Sommerfeld integrals) were adjusted to yield the fastest possible solution within an error 
of 5 X 10“^. Note that the last data points around /ci = Svr ~ 25.1 in Figure 9 (which is the last data 
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point presented for the layer-Green-function method) shows that, for such frequencies the WGF is 
approximately three orders of magnitude faster than the layer-Green-function method [18]. 



Figure 9: Gomputing times required by the WGF method (green line) and the layer-Green-function 
method [18] (blue line) to create the linear systems of equations resulting from the Nystrom dis¬ 
cretization of the relevant integral equations. 

The problem of scattering by the city-like structure depicted in Figure 10 is considered next. 
Figure 10 also displays the window function utilized in this example, which has been amplified 
by a factor 8 for visualization purposes. In contrast with the results presented previously in this 
paper, the case of TM-polarization is considered for this test. In order to properly account for the 
singular behavior of the fields near corners, the necessary graded meshes were generated utilizing 
the value p = 4 in the method described in [13]. Table 1 reports the computing times required 
to form the relevant system matrices for both the WGF method and the layer-Green-function 
method. Both solvers were optimized to produce a maximum error of 5 x 10“^ in the solutions of 
the integral equation, and the same computational grids were utilized to discretize the buildings 
for both methods. 

Table 1 compares the computing times required by the WGF method and the layer-Green- 
function method for two values of k 2 - In particular we note that, not only is the new method 
much faster than the previous approach, but also that the speed-up factor grows: a speed up 
factor in the hundreds for the value ^2 = 27r is doubled as k 2 is itself doubled to the value k 2 = dvr. 
Additionally, application of the layer-Green-function method in this context requires use of fictitious 
curves underneath each building [18] each one of which (curves) must itself be discretized, while 
the WGF method requires discretization of the ground between the buildings and in the region 
where the windowing takes place. In the present case the layer-Green-function method produced 
a system of 2384 unknowns while the WGF method produced a nearly identical sized system of 
2406 unknowns. At higher frequencies, the WGF method requires fewer unknowns than the layer- 
Green-function method, since, as demonstrated in Table 2, at higher frequencies the width of the 
windowing function can be decreased while maintaining accuracy. 

As an additional example we consider once again the city-like structure depicted in Figure 10 
but assuming an absorbing media in the ground and buildings: here we thus take ki = 27r and k 2 = 
47 r(l -|-i/100). Figure 11 demonstrates the convergence of both the naive windowing algorithm (6) 
and the full WGF method (10). The advantages provided by the full WGF approach can be 
appreciated clearly in this figure: in the naive method convergence near grazing is extremely slow 
while for the full WGF method the convergence is actually faster near grazing than for non-grazing 
configurations. In particular, the WGF method requires no more than 5 wavelengths of ground for 
a full four digits of accuracy, independently of the incidence angle. 
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Figure 10: City-like geometry and windowing function used. 



^2 

LGFM time 

WGFM time 

ratio 

TT 

27r 

588 s. 

3.07 s. 

192 

TT 

dvr 

3579 s. 

9.10 s. 

393 


Table 1: Computing times required by the layer-Green-function method and the WGF method to 
produce integral equation solutions with an accuracy better than 5 x 10“^ for the city-like geometry 
displayed in Figure 10. 


kx 

^2 

A 

TT 

27r 

6.5 

27r 

dvr 

3.5 

dvr 

Svr 

1.75 

87r 

IGvr 

1.1875 


Table 2: Extent of the windowed region required by the WGF method (10) to maintain an accuracy 
of 5 X 10“® in the approximation of the surface fields for the problem of scattering from a semi¬ 
circular bump of unit radius with various wavenumbers. The angle of incidence was taken to equal 
a = —7r/8 . 


For our last numerical example we consider an obstacle above the ground, but not connected 
to it, with a finite number of indentations under the ground level. Figure 12 displays the geom¬ 
etry under consideration, together with a selection of window function which yields an error of 
approximately 1% in the integral equation solution and corresponding near fields for a plane-wave 
illumination with incidence angle equal to a = —tt/S from the horizontal under TE polarization. 
Once again, as demonstrated in Figure 13 exponential convergence is observed as AjX grows. 
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Figure 11: Errors in the integral densities resulting from numerical solution by means of the layer- 
Green-function method (6) (left) and the WGF method (10) (right) for the city-like structure 
depicted in Figure 10, for various window sizes and angles of incidence—including extremely shallow 
incidences. Clearly, the WGF method computes integral densities with super-algebraically high 
accuracy uniformly for all angles of incidence. 



Figure 12: Scattering geometry containing a kite structure above a finite rectangular grating in 
an otherwise undisturbed planar ground. A windowing function large enough to produce an error 
smaller than 1% in the integral equation solution is shown along with the corresponding near fields; 
ki = 2 tt and k 2 = dvr. 


A Green function for a two-layer medium 

Consider the Helmholtz equation in the regions Hi = {(xi,X 2 ) G > 0} and CI 2 = {{xi-,X 2 ) G 

M^,X 2 < 0} with respective wavenumbers ki and k 2 - The Green function of the problem satishes: 

+ k]G = -6y 

^1x2=0+ ^\x2=0~ 

dG , _ dG , 

9x2 ^9x2 ^2=0 

and the Sommerfeld radiation condition at infinity, where 5y denotes the Dirac delta distribution 
supported at the point y. As is known G can be computed explicitly in terms of Sommerfeld 


in H 


•ji 


on {x 2 = 0}, 
on {x 2 = 0}, 


(16) 
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Figure 13: Errors in the integral densities resulting from numerical solution of (10) for the structure 
depicted in Figure 12 by means of the full WGF method, for various window sizes and angles of 
incidence—including extremely shallow incidences. Left: log-log scale. Right: semi-log scale. Once 
again we see that, the WGF method computes integral densities with super-algebraically high 
accuracy uniformly for all angles of incidence. 


integrals. To obtain such explicit expressions, given a fixed point y we define the functions ^j{x) = 
G{x,y), X G Qj. Expressing (pj as inverse Fourier transforms 


ipj{xi,X2) 


1 


(17) 


and replacing (17) in (16) a system of ordinary differential equations for the unknown functions pj 
is obtained which can be solved analytically. Two cases arise. For y G fli, the solution of the ODE 
system is given by 




$2{^,X2) 


g-7ih2-y2| / 1 _ jy\ g-7ih2+y2| 

271 + 271 

v{k2 — ki) 

(71 + z^ 72 )(l + I ') 71(71 + 72) ’ 

g-7102-3:2) / Ql2X2-liy2 g-7102-3:2) \ 

(1 -h 13)71 I 7i + ^12 (1 + ^ 3)71 j ’ 


(18) 


where -fj = — A:|. The determination of physically admissible branches of the functions 7j(0 = 


~ kj require selection of branch cuts for each one of the two associated square root 

functions. The relevant branches are — 37 r /2 < arg(^ — kj) < 7 r /2 for — kj and — 7 r /2 < 
arg(^ + kj) < 37 r /2 for + kj. Taking the inverse Fourier transform (17) of (pj and using the 
identity 


/ 


00 p-73 02-3/2! 


47r7j 


g3?(3:i-3/i) de = ^77i')(%|y-a^|). 


we obtain 


7^1 (®) =\Ho\ki\x - y|) + ^ (tT^) 

+ ^i{x,y), (19) 

H^ 2 {x) = - y\) + ^ 2 {x,y), 
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where the functions are given by 


^iix,y) = 


i>{k 2 — fcf) /■°° e '>'i(a; 2 +j/ 2 ) cos(^(xi — yi)) 


T^i'^ + v) Jo 


dC, 


1 ( g72a;2—71J/2 g7l(3:2—J/2) 

^2{x,y) =- 


TT 


71(72 + 71 ) (71 + ^ 12 ) 

cos(C(xi - yi))di, 


7i + ^12 (1 + i^)7i 

\ / 

Similarly, the solution of the ODE system for y G D 2 is given by 

j/g-7ia:2+72?;2 jy Q-'Y2{x2-y2) 


12 Q- 12 {x 2 -y 2 ) 

^l{C,X2) =— 77 --7- + 

(1 + U)j2 

Q-'r 2 \x 2 -y 2 \ 

$2{C,X2) = -7-h 


272 


7l + v ^2 (1 + 1^)72 

12 -l\ 6-721*2+^21 


12+1 


272 


+ 


I'ikf — k"^) e72(*2+j/2) 


(71 + I272)il + 12 ) 72(72 + 71 ) ' 

Taking inverse Fourier transform (17) we now obtain 

-^\k2\x - y\) + ^i{x,y), 

Hl^^\k2\x - y\) 


i 12 

+i(x) =7 


21 + 12 


+ 2 (x) =^-H^o\Mx - y\) + \ 


12-1 


12 + 1 


+ ^2ix,y), 


where the functions are given by 


^i(x,y) =- 


7T 


^ 00 / 072 ^ 2 - 71^2 Q-'l 2 (x 2 -y 2 ) 

) \ 7i + ^72 (1 + 12)72 


cos(7xi - yi))dC, 


U; n,\ _ ^(^1 ~ ^2) /■°° e^^(^^+^^)cos( 7 xi -1/1)) 

^ ’ 71(1 + 12 ) Jo 72 ( 71 + 72 ) (71 + 1 ^ 72 ) 


( 20 ) 


( 21 ) 


( 22 ) 


The gradient of the Green function is evaluated from the expressions above by differentiation under 
the integral sign. 
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